function PFnoX(H_t,Se,Phi1,y,n_agg,Mpf)

kap = 1

(T,n) = size(y)
n_s   = n-n_agg
s     = zeros(T,n_s)

# partition y
y_obs = y[:,1:n_agg]
s_obs = y[:,n_agg+1:end]

# partition Phi1
Phi_yy = Phi1[1:n_agg,1:n_agg]
Phi_ys = Phi1[1:n_agg,n_agg+1:end]
Phi_sy = Phi1[n_agg+1:end,1:n_agg]
Phi_ss = Phi1[n_agg+1:end,n_agg+1:end]

# partition Se
Se_yy = Se[1:n_agg,1:n_agg]
Se_ys = Se[1:n_agg,n_agg+1:end]
Se_sy = Se[n_agg+1:end,1:n_agg]
Se_ss = Se[n_agg+1:end,n_agg+1:end]

# Initialization
aParticles  = zeros(T,Mpf,n_s)
mWeights    = zeros(T,Mpf)
vLogLik     = zeros(T,1)

# Generate t=1 particles
vvs1 = [fMyMNDraw(s_obs[1,:],H_t[:,:,1]) for jj = 1:Mpf]
[aParticles[1,jj,:] = vvs1[jj] for jj = 1:Mpf]
mWeights[1,:] .= 1

# Forward Iterations
for t=2:T

    # println("Period t= $(t)")

    # Some expressions can be pre-computed because they don't depend on jj
    Pprime_ss_y_obs = Se_ss - Se_sy*inv(Se_yy)*Se_sy';
    F_yy            = Se_yy;
    F_ss_yy         = Pprime_ss_y_obs + H_t[:,:,t]
    kgain           = Pprime_ss_y_obs*inv(F_ss_yy);
    P_t_t           = Pprime_ss_y_obs - kgain*Pprime_ss_y_obs;
    ln_omega_t_fix  = -0.5*log(det(Pprime_ss_y_obs))+0.5*log(det(P_t_t));

    ln_omega_tilde_t = zeros(Mpf)
    aParticles_t     = zeros(Mpf,n_s)

    for jj = 1:Mpf

        # forecast state and yobs
        yprime = Phi_yy*y_obs[t-1,:]+Phi_ys*aParticles[t-1,jj,:];
        sprime = Phi_sy*y_obs[t-1,:]+Phi_ss*aParticles[t-1,jj,:];

        # condition distr of state on y_obs
        sprime_y_obs    = sprime + Se_sy*inv(Se_yy)*(y_obs[t,:]-yprime);
        # Pprime_ss_y_obs = Se_ss - Se_sy*inv(Se_yy)*Se_sy';

        # construct proposal density using KF updating step

            # forecast observation
            yobs_pred = yprime
            sobs_pred = sprime_y_obs;
            v_y       = y_obs[t,:] - yobs_pred;
            v_s       = s_obs[t,:] - sobs_pred;
            # F_yy      = Se_yy;
            # F_ss_yy   = Pprime_ss_y_obs + H_t[:,:,t]

            # updating step
            # kgain     = Pprime_ss_y_obs*inv(F_ss_yy);
            s_t_t     = sprime_y_obs + kgain*v_s;
            # P_t_t     = Pprime_ss_y_obs - kgain*Pprime_ss_y_obs;

        # draw particle values
        aParticles_t_jj = fMyMNDraw(s_t_t,kap*P_t_t)
        aParticles_t[jj,:]  = aParticles_t_jj
        # compute importance weights
        ln_omega_t_jj = (ln_omega_t_fix
            -0.5*(aParticles_t_jj - sprime_y_obs)'*inv(Pprime_ss_y_obs)*(aParticles_t_jj - sprime_y_obs)
            +0.5*(aParticles_t_jj - s_t_t)'*inv(P_t_t)*(aParticles_t_jj - s_t_t) )
        # compute incremental weights
        ln_omega_tilde_t[jj] = (ln_omega_t_jj -0.5*n_agg*log(2*pi) - 0.5*log(det(F_yy)) - 0.5*v_y'*inv(F_yy)*v_y
        -0.5*n_s*log(2*pi) - 0.5*log(det(H_t[:,:,t])) - 0.5*(s_obs[t,:]-aParticles_t_jj)'*inv(H_t[:,:,t])*(s_obs[t,:]-aParticles_t_jj) )

        # println("ln_omega_tilde_t(jj) = $(ln_omega_tilde_t[jj])")

    end

    # Likelihood increment:
    vHelp   = ln_omega_tilde_t .+ log.(mWeights[t-1,:]) #log of w-tilde * W
    c       = maximum(vHelp)
    vHelp   = vHelp .- c

    # println("vHelp(1) = $(vHelp[1])")

    vLogLik[t] = log(mean(exp.(vHelp))) + c

    # println("vLogLik(t) = $(vLogLik[t])")

    # Normalized mWeights
    vNormWeights = exp.(vHelp) ./ mean(exp.(vHelp))

    # Selection
    ESS                     = Mpf / mean( vNormWeights.^2 )
    doSelect                = ESS < Mpf/2

    if doSelect == 1

        #vSelectionIndices       = rand(Categorical(vNormWeights./Mbspf),Mbspf,1)
        vSelectionIndices   = resampleNYFED(vNormWeights)
        aParticles[t,:,:]   = aParticles_t[vSelectionIndices,:]
        mWeights[t,:]       .= 1

    elseif doSelect == 0

        mWeights[t,:]       = vNormWeights
        aParticles[t,:,:]  = aParticles_t

    end

    # s[t,:] = (mean*(vNormWeights.*aParticles[t,jj,:],dims=[1 2]))'

end

return s, vLogLik

end

#------------------------------------------------

function PFwithX(H_t,Se,Phi1,y,n_agg,densdraws,employdraws,K_exact,knots,xgrid,coef_lambda,coef_mean,Mpf)

kap = 1

(T,n) = size(y)
n_s   = n-n_agg
s     = zeros(T,n_s)

# partition y
y_obs = y[:,1:n_agg]
s_obs = y[:,n_agg+1:end]

# partition Phi1
Phi_yy = Phi1[1:n_agg,1:n_agg]
Phi_ys = Phi1[1:n_agg,n_agg+1:end]
Phi_sy = Phi1[n_agg+1:end,1:n_agg]
Phi_ss = Phi1[n_agg+1:end,n_agg+1:end]

# partition Se
Se_yy = Se[1:n_agg,1:n_agg]
Se_ys = Se[1:n_agg,n_agg+1:end]
Se_sy = Se[n_agg+1:end,1:n_agg]
Se_ss = Se[n_agg+1:end,n_agg+1:end]

# Initialization
aParticles  = zeros(T,Mpf,n_s)
mWeights    = zeros(T,Mpf)
vLogLik     = zeros(T,1)

# Generate t=1 particles
vvs1 = [fMyMNDraw(s_obs[1,:],H_t[:,:,1]) for jj = 1:Mpf]
[aParticles[1,jj,:] = vvs1[jj] for jj = 1:Mpf]
mWeights[1,:] .= 1

# Forward Iterations
for t=2:T

    # println("Period t= $(t)")

    # Some expressions can be pre-computed because they don't depend on jj
    Pprime_ss_y_obs = Se_ss - Se_sy*inv(Se_yy)*Se_sy';
    F_yy            = Se_yy;
    F_ss_yy         = Pprime_ss_y_obs + H_t[:,:,t]
    kgain           = Pprime_ss_y_obs*inv(F_ss_yy);
    P_t_t           = Pprime_ss_y_obs - kgain*Pprime_ss_y_obs;
    ln_omega_t_fix  = -0.5*log(det(Pprime_ss_y_obs))+0.5*log(det(P_t_t));

    # Prepare time-t cross-sectional data
    densdraws_t     = densdraws[1:N,t]./K_exact[t]
    employdraws_t   = employdraws[1:N,t]
    selecteddraws_t = densdraws_t[employdraws_t.==1]

    ln_omega_tilde_t = zeros(Mpf)
    aParticles_t     = zeros(Mpf,n_s)

    for jj = 1:Mpf

        # forecast state and yobs
        yprime = Phi_yy*y_obs[t-1,:]+Phi_ys*aParticles[t-1,jj,:];
        sprime = Phi_sy*y_obs[t-1,:]+Phi_ss*aParticles[t-1,jj,:];

        # condition distr of state on y_obs
        sprime_y_obs    = sprime + Se_sy*inv(Se_yy)*(y_obs[t,:]-yprime);
        # Pprime_ss_y_obs = Se_ss - Se_sy*inv(Se_yy)*Se_sy';

        # construct proposal density using KF updating step

            # forecast observation
            yobs_pred = yprime
            sobs_pred = sprime_y_obs;
            v_y       = y_obs[t,:] - yobs_pred;
            v_s       = s_obs[t,:] - sobs_pred;
            # F_yy      = Se_yy;
            # Pprime_ss_y_obs = Se_ss - Se_sy*inv(Se_yy)*Se_sy';
            # F_ss_yy   = Pprime_ss_y_obs + H_t[:,:,t]

            # updating step
            # kgain     = Pprime_ss_y_obs*inv(F_ss_yy);
            s_t_t     = sprime_y_obs + kgain*v_s;
            # P_t_t     = Pprime_ss_y_obs - kgain*Pprime_ss_y_obs;

        # draw particle values
        aParticles_t_jj = fMyMNDraw(s_t_t,kap*P_t_t)
        aParticles_t[jj,:]  = aParticles_t_jj
        # compute importance weights
        ln_omega_t_jj = (ln_omega_t_fix
            -0.5*(aParticles_t_jj - sprime_y_obs)'*inv(Pprime_ss_y_obs)*(aParticles_t_jj - sprime_y_obs)
            +0.5*(aParticles_t_jj - s_t_t)'*inv(P_t_t)*(aParticles_t_jj - s_t_t) )
        # reconstruct likelihood function for cross-sectional observations
        coef_t        = coefRecover(aParticles_t_jj', coef_lambda, coef_mean)
        ln_pdf_X      = -length(selecteddraws_t)*logspline_obj(selecteddraws_t, coef_t', knots, minimum(xgrid), maximum(xgrid))

        # compute incremental weights
        ln_omega_tilde_t[jj] = (ln_omega_t_jj -0.5*n_agg*log(2*pi) - 0.5*log(det(F_yy)) - 0.5*v_y'*inv(F_yy)*v_y
                             + ln_pdf_X)

        # println("ln_omega_tilde_t(jj) = $(ln_omega_tilde_t[jj])")



    end

    # Likelihood increment:
    vHelp   = ln_omega_tilde_t .+ log.(mWeights[t-1,:]) #log of w-tilde * W
    c       = maximum(vHelp)
    vHelp   = vHelp .- c

    # println("vHelp(1) = $(vHelp[1])")

    vLogLik[t] = log(mean(exp.(vHelp))) + c

    # println("vLogLik(t) = $(vLogLik[t])")

    # Normalized mWeights
    vNormWeights = exp.(vHelp) ./ mean(exp.(vHelp))

    # Selection
    ESS                     = Mpf / mean( vNormWeights.^2 )
    doSelect                = ESS < Mpf/2

    if doSelect == 1

        #vSelectionIndices       = rand(Categorical(vNormWeights./Mbspf),Mbspf,1)
        vSelectionIndices   = resampleNYFED(vNormWeights)
        aParticles[t,:,:]   = aParticles_t[vSelectionIndices,:]
        mWeights[t,:]       .= 1

    elseif doSelect == 0

        mWeights[t,:]       = vNormWeights
        aParticles[t,:,:]  = aParticles_t

    end

    # s[t,:] = (mean*(vNormWeights.*aParticles[t,jj,:],dims=[1 2]))'

end

return s, vLogLik

end
